Power spectra of self-organized critical sandpiles 



Lasse Laurson 1 , Mikko J. Alava 1 and Stefano Zapperi 2 3 

1 Helsinki University of Technology, Laboratory of Physics, HUT-02105 Finland 

2 CNR-INFM, SMC, Dipartimento di Fisica, Universita "La Sapienza", P.le A. Moro 2 00185 
Roma, Italy 

3 CNR, Istituto dei Sistemi Complessi, Roma, Italy 

Abstract. We analyze the power spectra of avalanches in two classes of self-organized 
critical sandpile models, the Bak-Tang-Wiesenfeld model and the Manna model. We show 
that these decay with a l// a power law, where the exponent value a is significantly smaller 
than 2 and equals the scaling exponent relating the avalanche size to its duration. We discuss 
the basic ingredients behind this result, such as the scaling of the average avalanche shape. 
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Sandpile models have been introduced almost twenty years ago as a paradigmatic 
example of self-organized criticality (SOC) (1), the tendency of slowly driven dissipative 
systems to display a scale free avalanche response. Such ideas have had an enormous impact 
in different fields, ranging from magnetic systems (|2|), superconductors (y|), mechanics (|4|;|5j), 
to geophysics and plasma physics including in particular the magnetosphere HQS). The 
influence also extends beyond physics, to for example biology (|9j), human (heart) physiology 
(Tioh - and cognitive processes or neuroscience (U). 

The reason behind this success lies in the wide variety of non-equilibrium systems 
displaying an avalanche response to an external driving. One of the primary aims of SOC 
was, originally, to explain the wide occurrence of 1/ f a noise in natural phenomena, through 
a direct relation between avalanche scaling and spectral properties (1). This idea was soon 
refuted when two groups 

a 

published works independently claiming that sandpile 
models should lead instead to a Lorentzian spectrum, that is decaying as 1/f 2 at large 
frequencies. The theoretical arguments were supported by numerical simulations on relatively 
small system sizes (|12t 
non-critical sandpiles ( 



1: 



). A non-trivial 1 / /"-decay in power spectra has only been found in 
, or by an alternative definition of the noise signal dl5l) , but not in 
standard cases such as in the original Bak-Tang-Wiesenfeld (BTW) model HI) and stochastic 
Manna model (fl6L 

Sandpile models represent a useful idealization of avalanche propagation, capturing 
the main ingredients behind this process: a slow external driving, a local threshold - or 
non-linearity - for the dynamics and a dissipation mechanism. While a complete exact 
solution of sandpile models is possible only in some particular cases , the origin of the 
scaling behavior is now well understood in the realm of non-equilibrium critical phenomena 



19c 120). Systems presenting a transition from an absorbing state to a moving phase, or 
similarly a depinning transition d21t l22l). can be turned into SOC under a suitable combination 
of a driving and a dissipation mechanism Jl8; 19; 20; 21). Conversely, criticality in sandpile 



models can be related to an underlying depinning critical point (|23t I24|) . The scaling of 
the power spectrum (PS) in sandpile models can be contrasted to avalanche induced crackling 
noise, which is typically characterized by a power law distribution of amplitudes and by a non- 
trivial 1/ f a spectrum (1251) . The most studied condensed-matter examples include Barkhausen 
noise in ferromagnets (2) and acoustic emission in fracture (0;EI) and plasticity Jif]). 

In this letter, we show that, notwithstanding previous beliefs, classical sandpile models 
display non-trivial 1/ f a spectra, a < 2 depends on the model and dimensions. We compute 
by numerical simulations the avalanche spectrum of two classes of sandpile models: the 
original two dimensional BTW sandpile model and the stochastic two-state Manna model, 
in one, two and three dimensions (Id, 2d, 3d). These two models are now known to be in 
different universality classes. A further difference between the two classes of models is that 
stochastic sandpiles obey finite size scaling while the BTW model displays multiscaling (E^f) . 
We find that the power spectrum decays as P(f) ~ f~ a , with a = 1.59 ± 0.05 for BTW and 
a = 1.44 ± 0.05, a = 1.77 ± 0.05, a = 1.9 ± 0.1 for Manna in d = 1, 2, 3, respectively. 

The central idea as to why SOC models can exhibit varying a, with the details depending 
on the dimension and universality class, is based on self-affine fractal dynamics. Consider 
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the time series V(t), which records the number of "topplings" (local relaxation events) taking 
place in the sandpile during each parallel update of the whole lattice, one such update defining 
the unit of time. An avalanche is defined here as a connected sequence of non-zero values of 
V(t). If the average size (i.e. the total number of topplings) of such avalanches of duration 
T scales as (s(T)) ~ T 7s * and the dynamics is self-similar, then the average avalanche shape 
should follow 

V(T,t) = T^- 1 f shape (t/T), (1) 

where f s hape( x ) is a scaling function d28l) . The form of the power spectrum is obtained 
averaging the energy spectrum E(f\s) of avalanches of size s, which in general scales as 

E(f\s) = s 2 g E {f^s), (2) 



where Qe{x) is another scaling function ( 13; 28). The power spectrum can then be written as 
P{f) = I D(s)E(f\s)ds, where D(s) ~ s~ T is the probability distribution of avalanche sizes 
and the integral is bounded by the upper cutoff s*, so that 

P(f) = r^ 3 - T) f ^ dxx 2 - T g E (x). (3) 



If the in teg ral in Eq. © is convergent, we obtain a = 7 S i(3 — r) (as originally derived 
in Ref. (29)). In the opposite case, the final result crucially depends on the asymptotic 
behavior of Qe{x). Kertesz and Kiss assumed (^(x) oc 1/(1 + x 2 ^ st ), obtaining a = 2 
Hil). Jensen et al. approximate the avalanche shape with a box function, which implies 
9e{x) oc (1 — cos(x 1//7st ))/x 2 / 7st , yielding again a = 2 I12L More recently, Kuntz and 
Sethna (28) noticed that if the toppling dynamics in the avalanche is a local process, the 
released energy is an extensive function of the size s, or E(f\s) ~ s. From Eq. © it thus 
follows that g E (x) ~ A/x. This implies that for r < 2 (which is the case for sandpile models), 
the integral in Eq. © is dominated by the frequency dependent upper cutoff, yielding a = ^ st . 

Here we analyze numerically the above set of four test-cases. We measure the shape 
of the pulse associated to an avalanche in the models, the scaling behavior of the avalanche 
size for a given duration and compute the power spectra. Sandpile models are defined on a 
d— dimensional hypercubic lattice. On each site i of the lattice the height is an integer variable 
Z{. At each step the system is driven, a grain is dropped on a randomly chosen site raising its 
height by one unit (z^ — > + 1). When one of the sites reaches or exceeds a threshold z c 
a "toppling" occurs: Zi = Zi — z c and Zj — Zj + 1, where j represents the nearest neighbor 
sites of site i. In the BTW model z c = 2d and each nearest neighbor receives a grain after the 
toppling of the site i. In the Manna model z c = 2 and therefore only two randomly chosen 
neighboring sites receive a grain. A toppling can induce nearest-neighbor sites to topple on 
their turn and so on, until all the lattice sites are below the critical threshold. This process 
defines an avalanche. We use parallel dynamics, meaning that every overcritical site topples 
when the lattice is updated, one such update defining the unit of time. The slow driving 
condition implies that grains are added only when all the sites are below the threshold. Grains 
can leave the system from the open boundaries. After a transient, sandpile models reach a 
steady state with avalanches of all sizes. Here we consider linear system sizes ranging from 
L = 1024 to L = 16384 in Id, from L = 64 to L = 2048 in 2d and L = 32 to L = 256 in 3d. 
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Power spectra are measured considering the signal produced by the number of toppling 
events V(t) as a function of time. Waiting times between avalanches have no effect on 
the scaling of the high frequency parts of the power spectra reflecting avalanche spectral 
properties (i.e. frequencies corresponding to time scales smaller than that of the duration 
of the longest avalanche). We have checked this e.g. by inserting a constant number of 
zeros between every successive pair of avalanches, and by means of a slow continuous drive, 
corresponding to a Poisson distribution of waiting times. Therefore in what follows, we 
simply join the avalanches one after the other. In Figs, [l] and El we display the power spectra 
P(f) and (s(l/T)) of the Manna and BTW models for different system sizes. While for 
the smallest lattices the power spectrum might be fitted by a Lorentzian, for larger system 
sizes the tails are definitely not scaling as 1/ f 2 . Instead, by fitting to the scaling parts of the 
power spectra (frequencies higher than the one corresponding to the duration of the longest 
avalanche and lower than inverse of a cross-over time after which the avalanches will have a 
self-similar structure), we find a = 1.59±0.05 for BTW and a = 1.44±0.05, a = 1.77±0.05 
and a — 1.9 ± 0.1 for Manna in Id, 2d and 3d, respectively. Note that the results are contrary 
to those in dl2c Il3l) . whose results were obscured by the small system sizes reachable at 
the epoch. Instead, at least for the Manna model, the scaling of the power spectra follows 
quite nicely that of (s(l/T)) ~ (1/T)~^, with = 1.44 ± 0.05, 7 st = 1.73 ± 0.05 and 
7 st = 1.9 ± 0.1 in Id, 2d and 3d, respectively Q. The PS of the BTW has a "bump" for 
small frequencies (Fig. 12), which shifts to still smaller ones with increasing L. Furthermore, 
for BTW, (s(l/T)) exhibits slight curvature even for large avalanches, but still the agreement 
is fair. 

In order to check that the observed results follow directly from the derivation outlined 
above, we compute the energy spectrum E(f\s). The results reported in Fig. |3] confirm 
the scaling behavior predicted by Eq. [2] with g E ~ 1/x. For the avalanche shape, Eq. (U]), 
the models show slightly different properties. Originally, Ref. (H2I) employed a simple box 
function to approximate f shape- This form is very far from the correct one, as it is shown in 
the inset Fig.|3j While the avalanche shape of the Manna model is symmetric, a more detailed 
look at the BTW model reveals, that its avalanche shape can not be rescaled as a function 
of duration T. The avalanches develop slowly an asymmetry, which could be related with 
the observations of multiscaling in the BTW model (27). For the stochastic Manna sandpiles, 
in Id the assumption of scale-invariance holds nicely while in the 2d and 3d cases relatively 
small avalanches show a cross-over behavior so that only around T ~ 100 the scaling regime 
is reached. This naturally implies corrections to scaling, but nevertheless a = 7,^ holds rather 
well (and in particular a < 2). 

The non-Lorenzian PS also persists if such systems are driven slowly as long as the 
avalanches do not overlap much, corresponding to timescales shorter than the inverse one 
related to the drive frequency. It is theoretically important and intriguing that such a relation 
can be established between the sandpile critical exponents and the PS one. In spite of the 
relation between SOC and non-equilibrium phase transitions, we still lack for analytical 
pred ictions for the critical exponents in d < 4, d = 4 being the upper critical dimension 
(20). Thus, for d > 4 mean-field exponents are valid, so that 7 st = a = 2. In two or three 



Power spectra of self-organized critical sandpiles 



5 



dimensions, however, we would generally expect that a < 2, and thus in many real physical 
systems the expectation would be the same. 

To summarize, self-organized criticality leads rather generally to power- spectra that 
exhibit 1 / /"-noise with a < 2. This calls for, perhaps, a re-evaluation of experimental results 
in many cases ranging from large systems met in solar and astro-physics to the laboratory 
and carrying over to the understanding of brain dynamics dill) , biology and so forth. In other 
words, a < 2 does not imply the absence of SOC, but instead may indicate exactly the 
contrary. There are also many theoretical issues that open up, including higher-order power 
spectra (13 ll) . In any case, we may safely conclude that the central point of the original paper 
by Bak, Tang and Wiesenfeld (i.e. the relation between avalanche scaling and non trivial 
power spectrum) was ultimately correct, contrary to what has been believed to be true from 
almost the beginning of the research on SOC and its applications to various phenomena. 
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f, 1/T 



Figure 1. The scaling of the power spectra (symbols connected by a line) and (s(l/T)) 
(symbols without a connecting line) of the Manna sandpile model for different system sizes, 
Manna Id (top), 2d (middle) and 3d (bottom). In all three cases, the slope of the spectrum is 
significantly different from 2 (dashed line). 
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f. 1/T 



Figure 2. The scaling of the power spectrum (symbols connected by a line) and (s(l/T)) 
(symbols without a connecting line) of the BTW model in 2d for different system sizes. The 
slope of the high frequency part of the power spectrum is again significantly different from 2 
(dashed line). 




Figure 3. Main figure: the energy spectrum for different sizes s is collapsed according to 
Eq. 13 and the scaling function decays as 1/x. Inset: the collapse of the average avalanche 
shape , avalanche durations ranging from T = 200 to T — 500. All data are for the 2d Manna 
model. 
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